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ABSTRACT 

We introduce a versatile class of prototype dynamical systems for the study of complex bifurcation cascades of limit cycles. 
Including bifurcations breaking spontaneously a symmetry of the system, period doubling bifurcations and transitions to chaos 
Induced by sequences of limit cycle bifurcations. The prototype system consist of a 24-dlmenslonal dynamical system with 
friction forces /(F(x)) functionally dependent exclusively on the mechanical potential V(x), which Is typically characterized, 
here, by a finite number of local minima. 

We present examples for d = 1,2 and simple polynomial friction forces f{V), where the zeros of f{V) regulate the relative 
Importance of energy uptake and dissipation respectively, serving as bifurcation parameters. Starting from simple Flopf- and 
homocllnic bifurcations, complex sequences of limit cycle bifurcation are observed when energy uptake gains progressively 
In Importance. 


Introduction 

The notion of a prototypical dynamical system is used whenever a certain system allows to study a certain relevant phe¬ 
nomenon, being at the same time simple enough that it is amendable for numerical and (at least partial) analytic investigations, 
without showing additional properties which could obfuscate the study of the prime phenomenon of interest.'"^ Additionally, 
their dynamical behavior can often be understood in terms of general concepts, such as energy balance, symmetry breaking, 
etc. 

Examples of prototype systems are the normal forms of standard bifurcation analysis^’ ^ and classical systems, like the 
Van der Pol oscillator"' or the Lorenz model,® which have been of central importance for the development dynamical system 
theory. As an example we consider the Lienard equation 

x + f{x)x + g{x)=Q, (1) 

a generic adaptive mechanical system, which includes the Van der Pol oscillator and the Takens-Bogdanov system.''-^ The 
periodically forced extended Lienard systems with a double-well potential have also been studied by many authors (see e.g. 
the double-well Duffing oscillator^”"). 

In this paper we propose a new class of autonomous Lienard-type systems, which allows to study cascades of limit cycle 
bifurcations using a bifurcation parameter controlling directly the balance between energy dissipation and uptake, and hence 
the underlying physical driving mechanism. Though there are a range of construction methods for dynamical systems in 
the literature (see e.g.'^”*^), they tend however to assume a somewhat abstract view, such as starting from implicitly defined 
manifolds, or involve mathematical tools accessible only to researchers with an in-depth math training. In contrast to these 
methods, we provide here a mechanistic design procedure, based on the construction of attractors by the use of potentials, 
a concept often used in many interdisciplinary fields (e.g. in modeling cardiovascular systems'® or for solving optimization 
problems'®), making it easily accessible and implementable for other scientific communities (such as neuroscience, biology 
etc.) as well. 

As an introductory example for the role of balance between energy uptake and dissipation in both local and global bifur¬ 
cations we reconsider the Bogdanov-Takens system. 


x={x-^)x-V'{x), 


x = y 

y={x-ll)y-V'{x) , 


( 2 ) 


which is often used as a prototype system for homoclinic bifurcations.^ Here, the mechanical potential is a third order 
polynomial, as illustrated in Lig. 1 . The friction force is directly proportional to the velocity y and fixpoints of (2) correspond 
hence to the minima and the maxima of the potential V (.r). 
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Figure 1 . The Bogdanov-Takens system (2) with the potential function V(x) =x^/3 — jl and a friction term x— jJ. (top 
row), and its generalization (4) to a friction term pi —V{x) (bottom row), compare Eq. (6). Left column: The potential 
function together with the color-coded regions of respectively energy dissipation and uptake, compare Eq. (3). Right column: 
The phase planes at the respective homoclinic bifurcation points, with the respective unstable foci and the saddles by open 
circles. The green and blue trajectories are the stable and unstable manifolds and the red trajectory corresponds to the 
homoclinic loop. 


The dynamics of the Bogdanov-Takens system is controlled by the parameter p, defining, in terms of the mechanical 
energy E, the regions of dissipation and energy uptake in the potential valley, 

E = {x-p)y\ E = ^-+V{x), (3) 

compare Eig. 1. The region x > p of energy uptake increases when the bifurcation parameter p is decreased, leading to two 
consecutive transitions. Initially the potential minimum becomes repelling, undergoing a supercritical Hopf bifurcation and a 
stable limit cycle emerges. Decreasing p further the extension of the limit cycle increases, merging at Pc with the stable and 
unstable manifolds of the saddle, resulting in a homoclinic bifurcation. 


Results 

The key mechanism leading to a bifurcation in the Bogdanov-Takens systems is the availability of a parameter allowing to 
change the balance between energy uptake and energy dissipation along limit cycles. Our aim is to generalize this idea to the 
case of mechanical systems characterized by an arbitrary number of potential minima. Eor this purpose we consider with 


x = y 

y = /(y(x))y-Vy(x) 


£ =/(y(x))y2 


(4) 


a 2(i—dimensional system, with d spatial coordinates and with friction forces depending via f(V (x)) functionally only on the 
mechanical potential V (x), allowing a fine-tuned control of the energy dissipation and uptake around the respective potential 
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Figure 2. Left column: Double well potential, as defined by Eq. (7) and with xi 2 = ±1, Ei 2 = 0, zi .2 = 1 and 2 = 1- The 
regions of energy dissipation E <0 and uptake £ > 0 are color coded. For the friction functions (6) we used fi, with a = 1 
{top row), f 2 , with 112 = 0.6 and a — 5 {middle row), and / 3 , with jj ,2 = 0.3, = 0.6 and a = 5 {bottom row). Right column: 

Bifurcation diagrams of the respective generalized Lienard systems (4), as a function of pi. All other Pi (when present) are 
kept constant. Stable/unstable hxpoints or limit cycles are denoted by continuous/dashed curves respectively. Black lines are 
fixpoint lines, while the maximal/minimal amplitude of v in a cycle is denoted with red/green color. H points denote 
Hopf-bifurcations, HO corresponds to homoclinic bifurcations of a saddle, SNC points denote saddle node bifurcations of 
limit cycles. The dotted, dashed and continuous vertical gray lines are just guides for the eyes. 


minima. A well known example of the system of type (4) is the Van der Pol oscillator 

x-e{l-x^)x+x = 0, f{V) = e(\-2V), = y ’ (5) 

for which the regions of energy uptake and dissipation remain hxed, with e regulating the overall influence of the velocity- 
dependent force. 
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Figure 3. Flow diagrams for the systems presented in Fig. 2, using respectively linear / quadratic / cubic friction functions 
f\(y) ! fi^y) ! h(y) itop/middle/bottom row). The values = 0.1/0.2, =0.1/0.11 and = 0.032/0.04 for 

the respective fii are indicated by arrows in the corresponding bifurcation diagrams in Fig. 2. 


The simplest generic class of friction functions f{V) entering (4) are polynomial; 

/i(y) = -a(F-Mi), f2iV) = -aiV-pi){V-p2), h{V) =-a{V - )ii){V - p2){V - (6) 

where a regulates the overall strength of the friction and where the individual pi < 112 < fla are the respective zeros, the points 
at which dissipation changes to anti-dissipation and vice-versa, compare Fig. 2. 

When using /j (V) and the mechanical potential V{x) = x^/3 — x^/2, the resulting flow in phase space is equivalent to the 
one of the Bogdanov-Takens system (2), as shown in Fig. 1 . 

Generalized mechanical potentials 

We are interested in using (4) as prototype dynamical systems, especially for the case of non-trivial mechanical potentials V (x) 
having an arbitrary number M of local minima. One could in principle consider higher-order polynomials for this purpose. 
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Figure 4. Stable limit cycles and chaotic orbits of the Lienard prototype system (4) in a two-dimensional symmetric double 
well potential V (x) (color coded, as defined by Eq. (9)) and with a linear friction term fi {V (x)) = 0.5 (/Zi — y(x)). The 
bifurcation parameter/Xi is 0.1,0.15,0.25,0.265,0.2698,0.3 from (a) to (/). In (d) the four limit cycles can be mapped into 
each other by using the symmetry operations ffi 2 or 03 4, as discussed in the Methods section. For (e) only a single of the 
four stable limit cycles is shown. This needs to circle the two potential minima twice in order to retrace itself. In (/) an 
example of a chaotic trajectory is given. 


however they do not allow to control the overall height of the potential and the relative width of the local minimal in as simple 
fashion. 

For this purpose we use throughout this study potential functions of the kind 

^(x) =n + g„(z) =tanh(zVz«^) , (7) 

n \ Pn / 

where the > 0 determine the half-width of the respective local minima and where the satisfy the self-consistent condition: 



since g(0) = 0. For deep minima, with {zi+Zj) <C |x; — x^ |, the positions and the heights of the local minima are close to x„ 
and V„ respectively. We found, that a relative accuracy of 10^^ for Vj, can already be achieved in general after three or four 
iterations. 

Limit cycie bifurcation cascades 

The system of type (4) allows to describe complex cascades of limit cycle bifurcations and in Fig. 2 we show some illustrative 
examples using a symmetric double-well potential and linear / quadratic / cubic friction functions fi{V) / / 2 (V) / fitiV) 
respectively, see Eq. ( 6 ). We used numerical methods'^ to obtain the respective full bifurcation diagrams, with solid/dashed 
lines denoting stable/unstable fixpoints and limit cycles. The corresponding flow in phase space is illustrated in Fig. 3. 

For negative jii values the two fixpoints (±1,0) are stable, for the case of /i (V) and (V), and stable limit cycles evolve 
via two supercritical Hopf-bifurcation. For fziV), on the other side, a sub-critical Hopf bifurcation is observed at /Ti = 0. 
The respective stable/unstable limit cycles merge for fxiV) and / 2 (V) in a homoclinic bifurcation, whereas a more complex 
bifurcation diagram emerges for f^iV). Saddle node bifurcation of limit cycles are present for both / 2 (V) and /a (V). 
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Chaos via period doubiing of iimit cycies 

We consider now a prototype system (4) with a two-dimensional symmetric potential field V (x), 

V(x) =g(x-Xi)^(x-X 2 ), ^(z) = tanh(4z^/9) , (9) 

as defined in (7), having two minima Xi 2 = ±(1, —1), and a linear friction term fi(V) = 0.5(fi —V). Both diagonals in 
the (xi,X 2 ) are symmetries of the system, as discussed in the Methods section. In Fig. 4 we present examples of stable 
limit cycles and of a chaotic trajectory, as projected to the (xi,X 2 ) plane. In Fig. 5 the corresponding bifurcation diagram 
is presented, which shows Hopf bifurcations (FI), homoclinic bifurcations (HO), branching of limit cycles via spontaneous 
symmetry breaking (SSB), period doubling of limit cycles (PD) and a transition to chaotic behavior; 

H At = 0 the two potential minima become unstable, just as for the one-dimensional spatial system presented in Fig. 2, 
resulting in two equivalent super-critical Hopf bifurcations. We note that, as a result of the symmetric potential function 
(9), a second branch of limit cycles is created by the two Hopf bifurcations (see the discussion in the Methods section 
and the Supplementary Information). However, since in the parameter region of interest these limit cycles are mostly 
unstable, we have not investigated them in detail. 

HO At ~ 0.143 the limit cycles merge, as in Fig. 4(a) — (b), in a homoclinic transition. The limit cycle stays, however, 
exactly on the diagonal xi -|-X2 = 0. 

SSB At the first branching point of limit cycles, RiO.171 the symmetry with respect to the diagonal (1,-1) is sponta¬ 

neously broken, as in Fig. 4(b) — (c), with the two limit cycles still being symmetric with respect to the (1,1) diagonal. 
The latter symmetry is broken at the second branching point Ri 0.260, as in Fig. 4(c) — (d), with four symmetry 
related limit cycles being stable. 

PD For larger values of the bifurcation parameter /ti a series of period-doubling of limit cycles is observed, with the first 
occuri'ing at ~ 0.268, as in Fig. 4(d) — (e). The next period-doubling transition occurs at r; 0.270, as 
shown in Fig. 5. 

For reference we note, that the saddle of the potential is located at V (0,0) = 0.505, viz at a substantially larger value. 

For Hi > ~ 0.2705 we observe seemingly chaotic trajectories, as illustrated in Fig. 4(f). Studying the transition to 

chaos is not the subject of the present investigation and we leave it to future work. We presume however, that the transition 
occurs via an accumulation of an infinite number of of period-doubling transitions of limit cycles, similar to the ones observed 
for the Lorenz system** and for the Rossler attractor. *^’^** 

Our prototype system (4) is not genetically dissipative. We have evaluated the average contraction rate a, as defined by 
(22) in the Methods section, and present the results in Fig. 5. Phase space contracts trivially along the attracting limit cycles, 
but also, on the average, in the chaotic region, where the average Lyapunov exponent A becomes positive. A is negative 
also for Hi < 0^ when only stable fixpoints are present, vanishing for intermediate values of Hi, when stable limit cycles are 
present. The later is due to the fact, see Fig. 7 and the cotTesponding Methods section, that two initially close trajectories will 
genetically flow to the same limit cycle with the relative distance becoming then constant. 

For larger values of Hi > 0.322 the chaotic region transients into a phase of intermittent chaos as illustrated in Fig. 6, in 
which an extended quasi-regular flow along the (—1,1) diagonal is interseeded by a roughly perpendicular bursting flow. This 
behavior is, to a certain extend, reminiscent to a scenario of intermittent chaos,^* in which a strange attractor is embedded 
in a higher-dimensional space with partly unstable directions. We have, however, not investigated the observed intermittent 
dynamics in detail. 

Discussion 

We have proposed and discussed a prototype dynamical system (4) in which the friction forces f{V ) depend functionally 
only on the mechanical potential V (x). We have shown, that complex cascade of limit cycle bifurcation can be obtained even 
for two dimensional phase spaces when the friction function f{V ) alternates between regions of energy uptake and dissipation. 

We have also introduced a generic class of potential functions (7) which allow to define, in a relative straightforward 
manner, mechanical potentials with an arbitrary number of local minima with varying depth. Any other potential could be 
however used. For example one could study the biquadratic version 

y(x) -)> (x-Xl)^(x-X 2 )^ (10) 
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Figure 5. Top row: The numerically obtained bifurcation diagram for a two dimensional prototype Lienard system with 
symmetric double-well potential and a linear friction force, as for Fig. 4. The second branches of limit cycles emerging from 
the two destabilized minima (see Supplementary Fig. SI) are not shown here. One observes Hopf and homoclinic 
bifurcations (H and HO), branching of limit cycles via spontaneous symmetry breaking and period doubling (SSB and PD), 
as well a transition to chaos for p > 0.2705. The red/green lines indicate the maximal/minimal xi-values of the respective 
limit cycles. The blue curve is the second xi-maxima after period doubling. The right diagram represents a zoom-in of the 
transition to a chaotic region, indicated by the shaded green and red areas. Only the first two period doubling bifurcations are 
shown. Bottom row: The average Lyapunov exponent A and the contraction rate (7, calculated as described in the Methods 
section, for the corresponding jii parameter intervals. For the left figure Api — 0.001 parameter stepsize was used. 

Increasing the resolution more and more periodic windows (with A = 0) become visible, as shown on the right plot, where 
Ajj .1 is decreased with a factor of ten. 


of the potential (9) used in our study of chaotic behavior with the prototype system (4). We did not study in detail the 
bifurcation diagram for the potential function (10), presume however, that it would be similar to the results presented in Fig. 4 
and 5, having the same underlying driving mechanism in terms of a linear friction function f{V) {fli — V). 

We have shown that a simple double-well prototype system with two spatial dimensions (and with a four-dimensional 
phase space) shows both symmetry induced bifurcations of limit cycles together with a period-doubling of limit cycle transition 
to chaotic behavior. 

As a future perspective we note, that by changing the depth of the minima, one could control the order in which the fixpoints 
are going be destabilized, which might lead to other interesting phenomena. Adding an extra (maybe slow) dynamics to the 
positions or heights of the minima, the metadynamics of the attractors^^ may also be considered. In this case the p„ parameters 
should be recalculated in each time-step using the self-consistent equations (8), which proved to be fast enough for practical 
purposes. 

Models, for which the equations of motion are derived from higher order principles, provide promising results for the un¬ 
derstanding many different phenomena, such as the optimization hardness of boolean satisfiability problems'® or the complex 
dynamics of biological neural networks.^^’^^ Generally, these methods involve the construction of a generating functional, 
such as the cost function or energy functional,with the dynamics of the system being defined by a gradient decent rule. 
When all equations are derived from the same generating functional, the system corresponds mathematically to a gradient 
system for which the asymptotic behavior is determined by the stable fixpoints (nodes). Thus they can not produce limit 


7/11 


































0 5000 10000 15000 20000 25000 30000 - 4-20 2 4 

t X2 


Figure 6. Left: Quasi-periodic windows with different time scales in the time-series plot of xi variable for jii = 0.34. Right: 
Phase space plot of the trajectory projected to the (xi ,X 2 ) plane for the same parameter. The magenta curve corresponds to 
the fp = [4 • 10^, 12 • 10^] quasi-periodic time interval, denoted with the same color in the xi (t) plot. 


cycles or any oscillatory behavior. To by-pass this problem, usually additional equations of motions are dehned, derived either 
from a second generating functional to produce objective function stress^®’^** or from other considerations. In our model, the 
system has an inherent inertia, which leads to damped oscillations around equilibria (minima of the potential) in the presence 
of dissipation. By creating regions of antidissipation, stable oscillatory dynamics and chaotic behavior can also be achieved. 
Considering nonsymmetric, higher dimensional potential functions we expect to hnd an even richer set of dynamical behaviors 
(see the Methods section), a scenario worth be investigated in the future. 

Finally, we note that the concepts of dynamical systems theory, such as attractors, slow points and bifurcations have been 
used recently to understand phenomena of surprisingly diverse fields. For examples of some relevant applications we mention 
here the modeling of birdsong,^' migraine^^ dynamics and the control mechanisms developed for the movements of humanoid 
robots.Hence, we believe that prototype systems which allow the construction of models with predefined attractors in an 
intuitive manner could offer a useful tool for understanding the behavior of interesting interdisciplinary problems. 


Methods 


Hopf bifurcations in the prototype system 

2-dimensional prototype systems 

The local maxima of the potential function, i.e. where V” (x*) < 0, which are saddle points, separate the phase plane into 
different attraction domains with their stable manifold. Local minima with V”{x*) > 0 become, on the other hand, repelling 
focuses as a result of an Andronov-Hopf bifurcation, when dissipation changes to antidissipation in their neighborhood, having 
a simple pair of purely imaginary eigenvalues X 12 = (x*). 


4-dimensional prototype systems 

Analogously, the 


q ={xuX2,yuy2)^ 


yu = o, 


dV{xi,X2) 


dx 


1.2 


= 0 


fixpoints of the 4-dimensional prototype systems (4) correspond to critical points of the V {xi,X 2 ) potential function. Classifica¬ 
tion of the local minima and saddle critical points with respect to their stability can be achieved by evaluating the eigenvalues 
of the Jacobian of the system in terms of the Hessian of the potential function: 


-/(q*) = 


/o 0 1 
000 

di c a 
\c d 2 0 

where we have defined with 


1 

0 

a j 


a = f{V{x\,xl)), 


c = 


H{x\,X*2) 



d^v 

^xidx2 ’ 


di.2 = — 


d^V 


(9x] 2 


( 11 ) 


( 12 ) 
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the friction term and the second order partial derivatives of the potential at the respective critical points. Dehning with 


7± = ^- {d, +d2)±^{d,+d2f- 4{did2 - c2)^ , 7 ±eM (13) 

we can express the eigenvalues of the Jacobian 7(q*) as 

■^1 .2.3,4 ~ 2 —47±^ . (14) 

For general potential functions the local minima, dehned by the A = det(//) =did 2 — c^ >0 and p = tr(//) = —{di+ d2) < 0 
conditions (or equivalently by 7 ± > 0), undergo a Hopf bifurcation, when the f{V) friction term changes sign, i.e.: 

a=fiy{x*i^x*2))=0 ^ Xi2,3A = ±iy/^- (15) 

However, saddles of the potential function, i.e. A = det(//) < 0, are saddle type fixpoints of the dynamical system, having 
always a Ai positive eigenvalue, as 7 + > 0 and j- < 0 . 

Here we note that in case of the potential function (9), due to the symmetries one gets a double pair of imaginary eigen¬ 
values, since d\= d 2 and c = 0, and thus Eq. ( 1 3) yields Y+ — Y-- This results in a second branch of limit cycle solutions, not 
investigated in this paper, emerging from the Hopf-point. 


2d-dimensional prototype systems 

For arbitrary dimensions d one can express the Jacobian in terms of block matrices; 




f Od Id\ 

\-Hd aid)' 


(16) 


where a = f{V) and where Od and Id are the (J-dimensional zero and identity matrices. Hd = (//,.;(x*)) = 

Hessian matrix of the V (q) potential, evaluated for the respective q* = (x*,y*) hxpoint. 

To determine the eigenvalues of the Jacobian one has to solve the equation: 


^ is the 


det(y — Xld) 


—Xld Id 

-Hd {a-X)Id 


deX{-X{a - X)Id + Hd) = 0, 


(17) 


where we used the properties of square block matrices. By introducing 7 = A (a — A) on hnds with 


d 

det(7 - Xld) = det(//rf - Yh) = Yi) = 0 

i=\ 


(18) 


that the 2d eigenvalues A of the Jacobian can be expressed in terms of the d eigenvalues 7 of the Hessian matrix and hence 

A* = i (^a ± ^/a^-AYi) . (19) 

Consequently, at the local minima of the potential, i.e. when 7 > 0, a Hopf-bifurcation occurs, with XX = when the 

friction term a = f{V ) changes sign. For general potential functions this might lead to the birth of higher dimensional tori or 
several branches of limit cycle bifurcations. 


Symmetries of the 4-dimensional system 

The results shown in Figs. 4, 5 and 6 are for a 4-dimensional prototype systems (4) with a linear friction force fi{V), as 
dehned in ( 6 ), and a mechanical potential V (x) given by (9). The minima V (xi 2 ) = 0 of the potential, viz. xi = (-fl, — 1) and 
X 2 = (— 1, -I-1) are connected by the symmetry operations 



/ 0 

±1 

0 

0 ^ 


/ 0 

±1 

0 

0 \ 


±1 

0 

0 

0 


±1 

0 

0 

0 

O ’!, 2 = 

0 

0 

1 

0 

J 0 ’ 3,4 = 

0 

0 

0 

±1 


^0 

0 

0 

V 


^0 

0 

±1 

0 / 


of the system. Thus, if (xi,X 2 , 7 i, 72 ) is a solutions, then 






An 


^xA 

x ’2 


2C2 


x" 

^ CT34 

2C2 

/l 

Ai 

5 

y'! 

yi 

Vi) 


Kyi) 


Kyi) 


Kyy 


( 20 ) 


( 21 ) 


are also solutions. 
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Figure 7. The logarithmic growth rate (In(Ar)) averaged for 100 random initial conditions as a function of time for three 
qualitatively different types of dynamics: spiraling into a fixpoint (jj.i = —0.05), limit cycle oscillations (jj.i — 0.2) and 
chaotic behavior (/Xi = 0.3). Brown lines correspond to the best linear regression. In the first and last case, the line is fitted 
only to the first part of the trajectory. The dashed line indicates that the distance of the point pairs has reached the maximal 
accuracy of the integrator. 

Lyapunov exponent and contraction rate 

The local Lyapunov exponent A is determined from the growth rate of the distance Ar(f) = Aroe'^', between point pairs with 
an initial displacement, which we have taken to be Aro = 10^^. Measuring the Lyapunov exponent was started after a transient 
of ttr = 1.5 • 10"*^. Considering 100 random initial conditions the average Lyapunov exponent A is then given by the slope of 
the initial linear part of the (In(Ar)) curve (as given by the brown lines in Fig. 7). 

The contraction rate a, is defined as the average of local contraction rates along a set of trajectories F for different initial 
conditions: 

a = (^^J^V-fdsy (22) 

where L = Jpt/i is the length of the trajectory and f is the flow, viz the right-hand side of the evolution equations (4). o is 
negative for dissipative systems, in which the phase space contracts.^’^ 
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